################################################################################
# Libraries and Import
################################################################################
rm(list=ls())

# Libraries
library(sandwich)
library(tidyverse)
library(stargazer)

# Import
load(file="data/dat_analysis_JP.RData")
jp <- vignette

# Set seed
set.seed(30)

# Control number of simulations in randomization procedures
sims = 10

################################################################################
# Clean covariates
################################################################################
# Japan
jp$age <- as.numeric(as.character(jp$age_group))
jp$international_sys <- as.numeric(as.character(jp$international_sys))
jp$perception_fp <- as.numeric(as.character(jp$perception_fp))
jp$critique <- as.numeric(as.character(jp$critique))
jp$Japan_citizen <- as.numeric(as.character(jp$Japan_citizen))
jp$territory <- as.numeric(as.character(jp$territory))

jp$nationalism <- 0
jp$nationalism[jp$vignette_group == 1| jp$vignette_group == 2] <- 1
jp$repression <- 0
jp$repression[jp$vignette_group == 2| jp$vignette_group == 4] <- 1

################################################################################
# Split into two datasets by treatment group
################################################################################
# Split dataset by treatment group - Japan
nat_jp = jp %>% filter(!is.na(nat)) %>% rename(z = nat)
econ_jp = jp %>% filter(!is.na(econ)) %>% rename(z = econ)

###############################################################################
# Estimation of heterogenous effects: China/Japan by Party
###############################################################################
# Add DPJ indicator
nat_jp$dpj <- with(nat_jp, ifelse((party == 6), 1, 0))
econ_jp$dpj <- with(econ_jp, ifelse((party == 6), 1, 0))

# Define ruling party in Japan as LDP only
nat_jp$party <- with(nat_jp, ifelse((party == 5), 1, 0))
econ_jp$party <- with(econ_jp, ifelse((party == 5), 1, 0))

# Create baseline tables: estimate models with no covariates
nat_jp_party = lm(D1 ~ z + party + z:party, data = nat_jp)
econ_jp_party = lm(D1 ~ z + party + z:party, data = econ_jp)

###############################################################################
# Effects by individual party (CDPJ vs. LDP)
###############################################################################
# DPJ vs LDP (Economic)
groupmeans_econ_dpj <- econ_jp %>%
  filter(dpj == 1) %>%
  group_by(z) %>%
  summarise(groupmean = mean(D1, na.rm = TRUE), size = n(), var = var(D1)) %>%
  mutate(party = "CDPJ")

groupmeans_econ_ldp <- econ_jp %>%
  filter(party == 1) %>%
  group_by(z) %>%
  summarise(groupmean = mean(D1, na.rm = TRUE), size = n(),var = var(D1)) %>%
  mutate(party = "LDP")

# Calculate standard errors
groupmeans_econ_dpj$se <- sqrt(groupmeans_econ_dpj$var/groupmeans_econ_dpj$size)
groupmeans_econ_ldp$se <- sqrt(groupmeans_econ_ldp$var/groupmeans_econ_ldp$size)

# Combine groupmeans estimates
groupmeans = rbind(groupmeans_econ_dpj, groupmeans_econ_ldp)

# Create chart
ate <- ggplot(groupmeans, aes(factor(z), groupmean, fill = factor(z))) +
  geom_col(position = "dodge", 
           color = "black", 
           width = 0.75) + 
  scale_fill_manual(values = c("steelblue2", "seagreen3"),
                    guide = FALSE) +
  geom_errorbar(aes(ymin = groupmean - 1.96*se, ymax = groupmean + 1.96*se),
                position = position_dodge(),
                color="black", 
                size=0.5,
                width = 0.75) +
  facet_wrap( ~ party, ncol = 2) + 
  geom_text(aes(label = round(groupmean, 2)), nudge_y = 0.5, size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Treatment") +
  ylab("Group means of government approval") + 
  scale_x_discrete("", labels = c("0" = "No interference", "1" = "interference")) +
  theme_classic() +
  scale_y_continuous(limits = c(0, 5)) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(legend.position = "bottom")

ggsave("figs/A5. ate_japan_ldp_dpj_econ.pdf", ate, height = 5, width = 6)

###############################################################################
# Tables of heterogenous effects: China/Japan by Party
###############################################################################
stargazer(nat_jp_party, econ_jp_party,
          se = starprep(nat_jp_party, econ_jp_party),
          dep.var.labels = "Government support",
          covariate.labels = c("Government interference",
                               "Ruling party affiliate",
                               "Interference x ruling party affiliate",
                               "Constant"),
          column.labels=c("Japan: Nationalist", "Japan: Economic"),
          title = "Conditional Average Treatment Effects by Party (excluding Komeito)",
          label = "tab:nat_econ_party_no_komei",
          align = TRUE,
          header = TRUE,
          omit.stat = c("f", "rsq", "ser"),
          out = "tables/TA10. vignette_party_no_komei.tex")

###############################################################################
# Charts of heterogenous effects: China/Japan by Party
###############################################################################
# Extract ATE and SEs from regression models
nat_jp_ldp = lm(D1 ~ z, data = nat_jp, subset = party == 1)
nat_jp_non_ldp = lm(D1 ~ z, data = nat_jp, subset = party == 0)
econ_jp_ldp = lm(D1 ~ z, data = econ_jp, subset = party == 1)
econ_jp_non_ldp = lm(D1 ~ z, data = econ_jp, subset = party == 0)

ate <- vector(mode = "numeric", length = 4)
se <- vector(mode = "numeric", length = 4)
party <- vector(mode = "character", length = 4)
country <- vector(mode = "character", length = 4)

ate[1] = summary(nat_jp_ldp)$coefficients[2]
ate[2] = summary(nat_jp_non_ldp)$coefficients[2]
ate[3] = summary(econ_jp_ldp)$coefficients[2]
ate[4] = summary(econ_jp_non_ldp)$coefficients[2]

se[1] <- sqrt(diag(vcovHC(nat_jp_ldp, type = "HC")))[2]
se[2] <- sqrt(diag(vcovHC(nat_jp_non_ldp, type = "HC")))[2]
se[3] <- sqrt(diag(vcovHC(econ_jp_ldp, type = "HC")))[2]
se[4] <- sqrt(diag(vcovHC(econ_jp_non_ldp, type = "HC")))[2]

party[1] = "LDP"
party[2] = "Other"
party[3] = "LDP"
party[4] = "Other"

country[1] = "Japan: Nationalist"
country[2] = "Japan: Nationalist"
country[3] = "Japan: Economic"
country[4] = "Japan: Economic"

party_ate = data.frame(ate, se, party, country)

# Create chart
ate <- ggplot(subset(party_ate, 
              country %in% c("Japan: Nationalist", "Japan: Economic")), 
  aes(party, ate, color = party)) +
  geom_point(size = 3, alpha = 0.9) + 
  scale_color_manual(values = c("steelblue2", "seagreen3"), guide = FALSE) +
  geom_errorbar(aes(ymin = ate - 1.96*se, ymax = ate + 1.96*se),
                position = position_dodge(),
                color="black", alpha = 0.3,
                size=0.5,
                width = 0.3) +
  geom_text(aes(label = round(ate, 2)), nudge_x = 0.25, 
            size = 3, color = "black") +
  facet_wrap( ~ country, ncol = 2) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Party Affiliation") +
  ylab("Mean change in government approval") + 
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(legend.position = "bottom")

ggsave("figs/A6. ate_japan_party_color_no_komei.pdf", 
       ate, height = 5, width = 6)